library(specr)
library(fixest)
library(readr)
library(broom.mixed)
library(texreg)
library(dplyr)
library(fixest)
library(patchwork)
library(ggplot2)

# Cross sectional linear FE models that use aggregate data (2000-2019)

source("funs_and_constants.R")

# load in data
epr_yearly_group_clusters_h1 <- read_csv("data/epro_data_h1.csv")
epro_year_group_clusters <- read_csv("data/epro_data_h2.csv")
epro_year_group_clusters_geo <- read_csv("data/epro_geo_data_h3.csv")

mean_h1 <- epr_yearly_group_clusters_h1 %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_orgs, resource_and_agriculture, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)),
    none_one_or_both_nats_agri_char = Mode(as.factor(none_one_or_both_nats_agri_char), na.rm = TRUE)
  )

mean_h2 <- epro_year_group_clusters %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_orgs, n_ed_religions, hhi_rel, religious_segments_bin, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)))

mean_h3 <- epro_year_group_clusters_geo %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_non_intersection_group_polygons, multiple_polygons_bin, spatial_hhi, n_orgs, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)))




#######
## H1
######


full_model_h1 <- feols(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                    cluster = ~countries_gwid, data = mean_h1 %>% as.data.frame())

full_model_h1_max <- feols(share_disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                        cluster = ~countries_gwid, data = mean_h1)

full_model_h1_n_resources <- feols(disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                                cluster = ~countries_gwid, data = mean_h1)

full_model_h1_n_resources_max <- feols(share_disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups |
                                         countries_gwid,
                                    cluster = ~countries_gwid, data = mean_h1)


screenreg(list(full_model_h1, full_model_h1_n_resources, full_model_h1_max, full_model_h1_n_resources_max),
          stars = stars)

# old table
# texreg(list(full_model_h1, full_model_h1_max, full_model_h1_n_resources, full_model_h1_n_resources_max),
#           stars = stars,
#           custom.coef.map = coef_map,
#           custom.model.names = c("Disagreement (1/0)", "Prop. Disagreement", "Disagreement (1/0)", "Prop. Disagreement"),
#           table = FALSE,
#           include.proj.stats = FALSE,
#           file = "tables/h1_crosssection_lin.tex")

texreg(list(full_model_h1, full_model_h1_n_resources, full_model_h1_max, full_model_h1_n_resources_max),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:2, "Prop. Disagreement" = 3:4),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h1_crosssection_lin.tex")



# refit the same model but with resources and agriculture as a cat var to improve the effect
# plot below
full_model_h1_graph <- feols(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                       cluster = ~countries_gwid,
                       data = mean_h1 %>% 
                                    mutate(resource_and_agriculture = if_else(resource_and_agriculture == 1, "Both", "None/Single")) %>%
                                    as.data.frame())

effect_plot_bin = marginaleffects::plot_predictions(full_model_h1_graph,
                                                    condition = "resource_and_agriculture",
) +
  theme_minimal() +
  labs(x = "Natural Resources\nand Agricultural Production",
       y = "Disagreement")

effect_plot_cat = marginaleffects::plot_predictions(full_model_h1_n_resources,
                                  condition = "none_one_or_both_nats_agri_char",
) +
  theme_minimal() +
  labs(x = "Combination of Natural Resources\nand Agricultural Production",
       y = ""
       )


effect_plots = (effect_plot_bin | effect_plot_cat)

ggsave(plot = effect_plots,
       filename = "plots/effect_plots_h1.pdf",
       dpi = 300, device = "pdf", width = 15, height = 9, unit = "cm")


########
### H2
########

full_model_h2 <- feols(disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                    cluster = ~countries_gwid, data = mean_h2)

# n rels as IV
full_model_h2_n_rels <- feols(disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                           cluster = ~countries_gwid, data = mean_h2)


full_model_h2_herf <- feols(disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                         cluster = ~countries_gwid, data = mean_h2)


# full_model_h2_no1clusters <- gfeolser(disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid
#                        +(1+ nightlight_total_pc_log + n_tek_groups |gwgroupid) + (1+ nightlight_total_pc_log + n_tek_groups |countries_gwid), cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(n_orgs > 1), family = "binomial")

full_model_h2_max <- feols(share_disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                        cluster = ~countries_gwid, data = mean_h2)

full_model_h2_max_n_rels <- feols(share_disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                               cluster = ~countries_gwid, data = mean_h2)

full_model_h2_max_herf <- feols(share_disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                             cluster = ~countries_gwid, data = mean_h2)

screenreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf, full_model_h2_max, full_model_h2_max_n_rels, full_model_h2_max_herf),
          stars = stars,
          custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6))

texreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf, full_model_h2_max, full_model_h2_max_n_rels, full_model_h2_max_herf),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h2_crosssection_lin.tex")

###########
### H3
##########

full_model_h3 <- feols(disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                    cluster = ~countries_gwid, data = mean_h3)

full_model_h3_max <- feols(share_disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                        cluster = ~countries_gwid, data = mean_h3)

full_model_h3_bin <- feols(disagreement ~ multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                        cluster = ~countries_gwid, data = mean_h3)

full_model_h3_max_bin <- feols(share_disagreement ~  multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                            cluster = ~countries_gwid, data = mean_h3)

full_model_h3_hhi <- feols(disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                        cluster = ~countries_gwid, data = mean_h3)

full_model_h3_max_hhi <- feols(share_disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                            cluster = ~countries_gwid, data = mean_h3)

full_model_h3_prop_no_outliers <- feols(disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                                     cluster = ~countries_gwid, data = mean_h3)

screenreg(list(full_model_h3, full_model_h3_max, full_model_h3_hhi, full_model_h3_max_hhi, full_model_h3_bin, full_model_h3_max_bin), stars = stars)


texreg(list(full_model_h3, full_model_h3_bin, full_model_h3_hhi, full_model_h3_max, full_model_h3_max_bin, full_model_h3_max_hhi),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h3_crosssection_lin.tex")


##########
# graph for null results
##########


effect_plot_set = marginaleffects::plot_predictions(full_model_h3,
                                                    condition = "n_non_intersection_group_polygons",
) +
  theme_minimal() +
  labs(x = "Number of Group Settlement Areas",
       y = "") +
  coord_cartesian(xlim = c(0, 40))

effect_plot_rel = marginaleffects::plot_predictions(full_model_h2,
                                                    condition = "religious_segments_bin",
) +
  theme_minimal() +
  labs(x = "Several Religious Segments (0/1)",
       y = "Disagreement"
  )


effect_plots = (effect_plot_rel | effect_plot_set)

ggsave(plot = effect_plots,
       filename = "plots/effect_plots_h2_h3.pdf",
       dpi = 300, device = "pdf", width = 15, height = 9, unit = "cm")







